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Abstract In radar systems, tracking targets in low signal-to-noise ratio (SNR) 
environments is a very important task. There are some algorithms designed 
for multitarget tracking. Their performances, however, are not satisfactory in 
low SNR environments. Track-before-detect (TBD) algorithms have been de- 
veloped as a class of improved methods for tracking in low SNR environments. 
However, multitarget TBD is still an open issue. In this paper, multitarget 
TBD measurements are modeled, and a highly efficient filter in the framework 
of finite set statistics (FISST) is designed. Then, the probability hypothe- 
sis density (PHD) filter is applied to multitarget TBD. Indeed, to solve the 
problem of the target and noise not being separated correctly when the SNR is 
low, a shrinkage-PHD filter is derived, and the optimal parameter for shrinkage 
operation is obtained by certain optimization procedures. Through simulation 
results, it is shown that our method can track targets with high accuracy by 
taking advantage of shrinkage operations. 

Keywords Multi-target tracking • Track-before-detect • PHD filter 



1 Introduction 

In order to extract target measurements, traditional tracking methods apply 
a detection threshold at every scan. The undesirable effect of detecting the 
sensor data, however, is that in restricting the data flow, it also throws away 
potentially useful information. For high signal-to-noise ratio (SNR) targets, 
this loss of information is of little concern [1]. For low SNR targets, this loss 
of information could be critical for a radar tracking system. Therefore, some 
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new algorithms using unthresholded data are more advantageous than the 
traditional methods in tracking low SNR targets. 

The concept of simultaneous detection and tracking using unthresholded 
data is known in literature as track-before-detect (TBD) approach [I]. TBD al- 
gorithms could improve the performance of a tracking system, which has been 
investigated for surveillance radar |2]- In [3] and [4], the advantage of TBD 
methods is discussed and many TBD methods are reviewed and compared. 
As a batch algorithm using the Hough transform [5], dynamic programming 
[B] or maximum likelihood estimation [7 , TBD could be implemented. These 
techniques operate on several data scans and in general require large compu- 
tational resources [I]. 

As an alternative, recursive TBD method is based on a recursive single- 
target Bayes filter [T]. An extension of the particle filter to multitarget TBD 
is given in [5], and an improved approach is given in [5]. In this algorithm, 
a modeling setup is applied to accommodate the varying number of targets. 
Then a multiple model Sequential Monte Carlo based TBD approach is used 
to solve the problem conditioned on the model, i.e., the number of targets [TU] . 
This approach has proven to be very efficient in both single and multitarget 
case [3], though it restricts itself to the case in which the maximum possible 
number of targets is limited and known. 

Another extension of the single-target Bayes filter to multitarget TBD is 
based on a multitarget Bayes filter. Because a single-target Bayes filter is op- 
timal for a single target, to solve the problems introduced by multiple targets, 
the multitarget Bayes filter is proposed in [11]. In a multitarget Bayes filter, 
multitarget states and observations are modeled as random finite sets (RFS). 
This approach is a theoretically optimal approach to multitarget tracking in 
the framework of finite set statistics (FISST) [T^]. However, the multitarget 
Bayes filter has no practical utility without an approximation strategy. To solve 
this problem, the probability hypothesis density (PHD) filter [13] is proposed 
as a tractable and calculation-simple alternative to the multitarget Bayes filter 

The PHD is the first moment of the multitarget posterior probability. Un- 
der some assumptions (e.g. Possion Assumption), the PHD is an approxima- 
tion of the multitarget posterior probability. Therefore, the PHD filter can be 
an approximation of the multitarget Bayes filter. Although a PHD algorithm 
for TBD is proposed in jlOj , this approach ignores TBD measurements should 
be modeled by RFSs. In [TJ], multitarget TBD from image observations is 
formulated in a Bayesian framework by modeling the collection of states as a 
Multi-Bernoulli RFS. This work use the multi-Bernoulli update to develop a 
high precision multi-object filtering algorithm for image observations, although 
its adaptability of low SNR environment is needed to discussed in more detail. 

In our previous work [16], we use the RFSs to model multitarget TBD 
measurements and the collection of states. In this way, a traditional PHD 
filter [m and [T^ could be applied to multitarget TBD. Even though PHD 
could be an approximation of the multitarget posterior probability, the accu- 
racy of the algorithm is limited by some reasons, which are indicated in this 
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paper. Firstly, when the SNR is too low, the PHD cannot be a sufficient ap- 
proximation of the multitarget posterior probability because the fundamental 
assumptions are challenged. Furthermore, for multitarget TBD, the measure- 
ments of target and noise can hardly be separated, while the PHD filter is 
heavily dependent on the measurements of targets jl7l . In traditional tracking 
systems, to solve these problems, cardinalized PHD filter (CPHDF) is pro- 
posed in [TB]. However, CPHDF is inefficient in multitarget TBD, because 
the computational complexity of CPHDF is O (m'^) , where is the number of 
elements of the measurement set and is quite large for the TBD problem. 

In this paper, for extending a traditional PHD filter to be suitable for mul- 
titarget TBD, viewed from a different perspective, TBD can be regarded as a 
kind of classification problem of target and noise measurements. To enhance 
the difference between the target and noise measurements and to pursue better 
classification performance, the measurements need to be denoiscd. The thresh- 
old shrinkage algorithm |19] is an important method for image denoising. In 
general, the key points of threshold shrinkage are the following: the method 
of shrinkage (e.g., soft-threshold method) and the selection criterion of the 
threshold. In this paper, a shrinkage operation is adopted that is similar to 
the threshold shrinkage algorithm. The optimal parameter for the shrinkage 
operation can be obtained via certain optimization procedures. 

Furthermore, in this paper, some problems of multitarget TBD in parti- 
cle use are also discussed. The assumption of known SNRs is used frequently 
in traditional TBD algorithms 0, [8] and [9]. Multiple targets with similar 
SNRs are common assumption in the simulations of PHD algorithms using 
the amplitude feature, as in |15| and |20] . However, in practical use, multi- 
ple targets with different or unknown SNRs are common. In this paper, for 
the SMC implement of the PHD filter adopted, this problem could be solved 
by augmenting the SNR into target state, varying the method of generating 
predicted particles and adjusting the update operator. 

In this paper, the measurement of targets is modeled by a 'nail-like' model 
on range-Doppler maps because of the assumptions of the classical PHD filter 
in the framework of FISST. Recently, a classical PHD filter has been mod- 
ified to solve the problems of extended and group targets in [21] and [22], 
respectively. Therefore, the TBD measurement model of extended targets in 
the framework of FISST will be discussed in future work. 

The rest of this paper is organized as follows. In Section 2, the multitarget 
TBD problem is modeled by RFSs. In Section 3, the limitation of traditional 
PHD filter extension to TBD is investigated. In Section 4, a shrinkage option 
for the PHD filter is proposed with an optimal parameter. Tracking multiple 
targets with different or unknown SNRs is discussed in Section 5. Simulation 
results for the tracking system are presented in Section 6. Finally, we conclude 
the paper in Section 7. 
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2 Multitarget TBD RFS model 

2.1 State of the RFS model 

The target state is = [a;^, ifc, j/fe, yfc]^, where {xk,yk) and {ik,yk) are the 
position and velocity. Since there is no ordering on the respective collections 
of all target states, they can be naturally represented as a finite set. Xk is 
the multitarget state-set at time step fc, i.e., the set of unknown target states 
(which are also of unknown number). 

In a single-target scenario, the state is modeled as random vectors. 
Then, the multitarget state, including target motion, birth, spawning, can be 
described by RFS. For target motion, given multitarget state set Xk-i, each 
Xfc_i £ Xk-i either survives at time step k with probability ek\k-i {^k~i), 
and its transition probability density from Xk-i to Xk is fk\k-i{'^k \^k-i) ■ 
Therefore, the target motion is modeled as the RFS Sk\k-i{^k-i)- In the 
same way, when the RFS of target birth at time k is modeled by 7^, and the 
RFS of targets spawning from a target with Xfc_i is modeled by Bk\k-i (xfe-i) 
, the multitarget stat Xk is given by 



Xk 



u 



2.2 Measurement of the RFS model 

The measurements are measurements of reflected power as in [8] . Uk is white 
complex Gaussian noise with variance ctq. In this section, assume that the 
intensity of all targets is Ik , the SNR for the targets is defined by 

SNR = 10 \og{-J^)dB (2) 
2ctq 

The method to deal with multiple targets with different or unknown SNRs 
is discussed in section 5. 

Because of the assumptions made in the modeling process and the essen- 
tial difference between an RFS and a random vector, the two conditions that 
should be satisfied to model TBD measurements by a RFS are the following: 
there is no target that generates more than one measurement vector, and no 
measurement vector is generated by more than one target PT]. In summary, 
the measurement of targets should be modeled by a 'nail-like' model on the 
range-Doppler-Bearing maps (the ijk cell is defined by coordinat (r.i, dj ,bi)) 
as follow: 

R, D and B are the size of a range, the Doppler and the bearing cell, — 

Jixkf + iVkf, dk = j;'\tlf\. and bk - arctan (|^) . i?^, = [r, - f , + f ) , 
Di ^ [d, f + f ) and Bi = [k - f ,6; + f ). 
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Then the measurement model is 

|2 



f ^\h'^\^k)+nA (3) 



\ else 

This means the measurement of a target is like a nail on the range-Doppler- 
bearing maps. This 'nail-like' model is similar to the point target model in 
[5] and [6]. At step k, the measurement provided by the sensor consists of 
N = Nr X Nd X Nb measurements z^i^\ where Nr, Nd and Nt, are the number 
of range, Doppler and bearing cells. 

The number of z^j^^ is a constant N. However, the number of element of 
an RFS should be random and Poisson distributed for a PHD filter. Because 
weak signal information should be preserved by TBD algorithms, we should 
make sure all measurements produced by targets are included in the RFS. 
Therefore, a threshold could be chosen as follows: 

Ok = argmin/+°°go(4^')rf4^' 

Ok (5) 

s.t. poi^k) - P{zf > Ok) = /,t°°.9fe(4'' \^k)dzt > 0.99 
where the likelihood function for the target is 

l^fe ) = ^ {- )^^\ — I (6) 



and that for noise is 



Po(4^') = ^exp|-^.f|. (7) 

where Iq () is the zero order Bessel function and z'^^^ is assumed positive. The 
Poi'^k) is the probability detection. Because this algorithm is for TBD, it 
should be ensured that pD{x.k) ~ 1- 

Zk is the observation-set consisting of all measurements collected by all 
sensors at time-step k, no matter the measurement from the targets or from the 
noise. If z^-*' > 6k, let Zk — z^^^ . Then, the measurement RFS Zk is constructed 
by a subset of the z^"*' using a thresholding mechanism. The threshold can 
insure that p£)(xfc) sa 1 by ([S]). The threshold is a function of the SNR of the 
targets. The element in the RFS Zk is as follows: 

'Ll^[ri,dj,buZkf (8) 

and the likelihood function according 

fffc(Zfc |xfe) |xfc) (9) 



6 



pSK) = ^exp|-f^^ (10) 



At the same time, the sensor could be modeled by 



,^leZk (11) 



where the measurements of targets are modeled by RFS 0k{xk) and the noise 
is modeled as RFS K^- The elements in Ok{xk) are the measurements pro- 
duced by the targets. The elements in Kk are the measurements which are 
produced by the noise and bigger than Ok as well. It is shown that when the 
measurements of TBD are modeled by RFSs, multitarget TBD can be regarded 
as a kind of classification problem of target and noise measurements for one 
scan. This classification problem will be further analyzed in Section 4. 



3 The traditional PHD filter extension to TBD 

After multitarget TBD measurements and the collection of states are modeled 
by RFS, a traditional PHD filter is applied to multitarget TBD. This algorithm 
is reviewed in 3.1. However, the accuracy of the algorithm is limited by some 
reasons when the SNR is low, which is discussed in Section 3.2. 



3.1 The algorithm 

As proposed in [T3], the PHD prediction and update equations are presented 
as follows: 



i'felfe-l (x|^l:fe_l) = /efc|fc-l(C)/fe(x|C)i>fc-l|fc-l(Cl^l:fc-l)f^C+ 

J bk\k-i i^\C)Dk-i\k-i (C \Zi:k-i )rfC+7fc(x) ^ ^ 



Dk\k{^\Zi:k) = [1-Pd{-x.)] Dk\k-iW+ 



(13) 

where Dk\k (x {Zi-k ) the PHD is the density whose integral Jg Dk\k (x \ Zi-k )dx 
on any region S of state space is %.(S') = J \X D S\pk{Xk \Zi.,k )SX. bk\k^i (xfe |xfe_i ) 
and 7fe(xfe) denotes the intensity ol Bk\k-i (xfc_i) and Jl- at time k, and Kk{zk) 
is the intensity of Kk- Zi-k is the time-sequence of observation-sets. 

Because the PHD propagation equations involve multiple integrals that 
have no computationally tractable closed form expressions. Sequential Monte 
Carlo (SMC) methods are used to approximate the PHD in [12 . Let the update 
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PHD at fc — 1 step Z)fc_i|fe_i (x|Zi:fc_i) be represented by a set of particles 

R\.xii};:;\as 



p=i 

The predicted particles are generated by 



(14) 



qk{» x^^'ij p = i,...,L 



k-l 



''''' ^ l^^fc(») P = Lk-i + l,...,Lk-i + Jk 



(15) 



where qk{^ x^^^ ) and Vk{^) are proposal density. The predicted density is 



D 



k\k-l 



where 



,(p) 

^k\k-l 



p=i 


-i) 


(16) 








p = 1,... 


Lk-1 


(17) 



Set Pd{x) = 1, the update density will be 

-Dfe|fe(x|Zi:fc) 



A: - 1 + -^fc 

p=i 



where 



(18) 



E 



p=i 



,(P) ,(P) 



(19) 



According to the standard treatment of particle filter, resample |w^*'^''/nj., : 
toget{4-V%,xi-)}^';^ 



'/c|/c- 
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3.2 The limitation 

Because Kk{zi^) can be presented by the number of noise samples Xk multiply- 
ing the clutter probability density, the update operator (jl9p turns into 



9k (^fc ^fcji-i) 



^ Lk-i + Jk 



" ^kP*o (zk : (To) + E gk* U*k ) 
p=i ^ ' ' 



PS(-fc:-)-^cxp|-^^| (21) 

Equation (PT|) is the likelihood function for the measurement of noise in the 
measurements set, which denotes the clutter probability density in theoreti- 
cal filtering ([T^] and [13]), and the optimal choice of a should be the vari- 
ance of noise , as shown in (j20l) . However, for TBD applications, the SNR 
is extremely low. Hence, this setting of leads to the situation in which mea- 
surements generated by the target and noise are not separated correctly and 
heavy degradation of tracking performance is seen. In fact, the derivation of 
(|20p depends on the hypothesis that the number of noise sample is Poisson 
distributed. This means the clutter follows a Poisson distribution and is inde- 
pendent of target-originated measurements. In this way, the PHD can be the 
"best fit" approximations of the multitarget posterior probability [13]. The 
average number of noise samples included in RFS is 



Afe=iVexp<j-^^. (22) 



When TV is large and is relatively small, the Poisson distribution could 
be approximated by the distribution of a number of noise samples. However, 
this approximation will be destroyed when become sufficiently large, which is 
the case in low SNR scenarios, as in Table 1. 



Table 1 Aversus SNR (N = 2000) 
SNR(dB) 6 7 8 9 10 

Afc 1340 1098 707 344 143 



When (To in ([20)) cannot make PHD a sufficient approximation of the multi- 
target posterior probability, choosing the "optimal" becomes the key problem 
for extending the PHD filter to be applied in TBD. As discussed in the fol- 
lowing section, this selection can be considered as a shrinkage operation, and 
the optimal parameter for the shrinkage operation can be obtained via certain 
optimization procedures. 



9 



4 Shrinkage operation for TBD 

To determine the "optimal" and make tlie noise and targets distinguishable, a 
different perspective should be taken. As referred to in Section 2.2, TBD can 
be regarded as a kind of classification problem for target and noise measure- 
ments. In other words, when a measurement is obtained, should it be classified 
into a target class or a noise class? To solve this problem, the Fisher class sep- 
arability criterion, as a supplement of to the traditional Bayesian framework, 
was adopted in our design for the extension of the PHD filter. Furthermore, to 
enhance the difference between the target and noise measurements and pursue 
the better performance of classification, some attempts to reduce the noise are 
needed. The threshold shrinkage algorithm [H] is an important method for 
the image denoising. In general, the cores of threshold shrinkage are as fol- 
lows: the method of shrinkage and the selection criterion of the threshold. A 
shrinkage operation is adopted in this section that is similar to the threshold 
shrinkage algorithm to some extent. The optimal parameter for the shrinkage 
operation, which acts as the threshold, can be obtained via certain optimiza- 
tion procedures. The PHD filter with the shrinkage operation is called the 
Shrinkage-PHD filter. 



4.1 Fisher class separability criterion 

It is known that the Bayesian classifier is the optimal classifier when the pos- 
terior probability can be calculated. However, as indicated in Section 3.2, the 
PHD cannot be a sufficient approximation of multitarget posterior probability 
when the SNR of the targets is low. Therefore, another classification criterion, 
the Fisher class separability criterion, is introduced into our methods. 

The separation between the two classes of objects defined by Fisher is the 
ratio of the distance between the centers of classes to the scattering of the 
classes 1^: 



^between _ (M1 A^o) 
f]'^ ^2 I J*2 ■ 

"'within "l "T "O 



(23) 



where Hi represents the centers of classes and denotes the scattering 
of the classes. In our analysis, they are set to the mean and variance of the 
likelihood functions of the target and the noise, respectively. 

For targets, the mean and variance are defined as follows: 

/•oo />oo 

= / zugk{zk |xfc )dzk ~ I Zkgk{zk \^k )dzk = 2al + ll (24) 



j zlgk{zk\xk)dzk ~ / zlgk{zk\xk)dzk ^ ialial + ll) (25) 
Je Jo 

For the noise, the mean and variance are given by 
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= Mo + ^fc = + Ok (26) 

df =4^ (27) 
Therefore, the Fisher class separabihty is 




Simulations indicated that the correlation between S and the SNR of the 
targets is a positive, approximately linear relationship (as shown in Fig[T]). 
Therefore, when the SNR is low, the Fisher class separability is so small that 
the targets and false alarms could not be distinguished. In other words, we 
should enlarge the Fisher class separability of the target class and the noise 
class when the SNR of the targets is low. 



3.5 




[I , , , , , , J 

6 7 8 9 10 11 12 13 

SNR(dB) 

Fig. 1 Fig.l The Fisher class separability versus the SNR 



4.2 Shrinkage operation 

To enlarge the Fisher class separability between targets and noise, some kind 
of " shrinkage" operation is applied to the likelihood function of the noise. That 
is, the parameter CTq in and (P?]) is adjusted. If as < ctq is chosen instead 
of (jg itself, then Fisher class separability becomes 



{2al + Il-e-2air 
^al{al+ll) + Aat 
Obviously, Sg is larger than S whenever (Tj < ctq ■ The following question 
naturally arises: how does one choose the "optimal" as , which is the key 
parameter for the shrinkage process? 
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4.3 The optimal parameter for the shrinkage operation 

To demonstrate how influences the classification problem of target and noise 
measurements, we can define the Mahalanobis distance [M] of measurement 
Zk using the center of the noise class and the target class, respectively: 

c f ^ (zfc-6lfc-2gg)^ 

Sso{zk,as)^- — (30) 



4cr4 



and 



Ssiizk) ^ ^^T^r^TTTT- (31) 



Because the selection of the threshold is a difficult task for the threshold 
shrinkage algorithm [TH] in the field of image denoising, the selection of the 
optimal parameter for the shrinkage-PHD filter involves similar problems. For 
the threshold shrinkage algorithm, an oversized threshold may result in the loss 
of signal. In contrast, too much noise may remain if an undersized threshold 
is used. In the same way, for the shrinkage-PHD filter, preserving the weak 
target information and reducing the noise by tracking are the major challenges 
to be addressed for TBD problems. Therefore, when considering the two sides, 
the selection criterion should be designed as follows: 

af'^ — max <Ts 
s-^- I Jo' 9k{zk |xfe )dzk = /3 ^^'^^ 

y Ssq{Zs,(Ts) > Ssl{Zs) 

The first constraint represents the " shrinkage" . To preserve the weak target 
information and classify as many low SNR targets to the target class as pos- 
sible, the Mahalanobis distance of measurement with the center of noise class 
should be enhanced; hence, the variance of the noise class should be reduced. 

The second constraint refers to the fixed probability of target loss. The 
reason for this term is that for TBD problems, false alarms can be reduced 
by integration over time. Meanwhile, if target loss occurred, the integration of 
the information of targets is broken off. Hence, the cost of target loss is much 
larger than the cost of false alarms, which is the main difference between TBD 
and common radar detection. Moreover, for the PHD filter, a missed detection 
can result in loss of the track [T7]. Therefore, target loss should be fixed as /?. 
When the (3 is set, a critical value Zs is determined. 

The third constraint means that the measurements Zk determined by the 
second constraint should be classified as targets. When the SNR is low, Sso{zs,crs) 
always intersects Ssi (zs) , so the optimal solution u*^ to ((32|) is deduced, as in- 
dicated in Figl^l^a). With the increasing SNR, Ssoizg, (Js) becomes much larger 
than Ssi{zs) and thus we choose cr^ = ctq, as shown in Figl^^b). 

Note that although the shrinkage operation could highlight the target mea- 
surement and reduce the possibility of losing the target during tracking, the 
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noise might be incorrectly increased unavoidably. To minimize the influence 
of shrinkage on the classification of the noise measurement, the maximum 
cTg , which is consistent with the above constraints, is obtained. The optimal 
solutions to under different SNRs are listed in Table 2: 



Table 2 a^' (/3 < 0.05, (tq = 0.25) 



SNR{dB) 


6 


7 


8 


9 




0.48O-0 


0.6(To 


0.68(70 


0.72(70 


SNR(dB) 


10 


11 


12 


13 




0.76(70 


0.88(70 


0.96(70 


(70 
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The update operator of the shrinkage-PHD filter can be summarized as 
follows: 



, ,*(p) _ 



E 



(P) 

k\k-l 



-l^^- Xp*{zk : af ) + E gli 



Lk-i + Jk 



(33) 



Ap) 

^k\k^ 



(p) 

k\k-l 



where should be set according to Table 2. 



5 Multitarget TBD in practical use 



This section is to solve of the problem of tracking multiple targets with different 
or unknown SNRs. 

As proposed in Section 4, the optimal parameter for the shrinkage operation 
is closely related to the SNRs of the target, and the SNRs of the target are 
assumed to be known and similar. However, in practical use, multiple targets 
with different or unknown SNRs are common. The PHD is the first moment 
of the multitarget posterior probability, but not the posterior probability of a 
certain target. Therefore, the Ik in ^ and the cr*^ in Tableware difficult to 
determine. Because the SMC implement of the PHD filter is adopted, this 
problem could be solved by adjusting the method of generating predicted 
particles and the update operator. 

To describe how to deal with the different SNRs of multiple targets, we 
augment Ik into : 



Xfe = [Xfe,/fe] 



(34) 



where x^ — [xk ,Xk,yk,yk]- The particle of the target state is x 
^ip) 



(p) 



The X uj,_i is generated as ([15 



-(p) Ap) 



^(P) 
-^fe|fe- 



Qki* x[^2i) P= 1,. 

Vk{») P = Lk-1 + 1, 



, Lk~l 
■ , Lk~l 



Jk 



(35) 



The /i^-* is generated by the assumption of a uniform distribution in the 



range of 1^ , where the SNRs of all of the targets are 
and corresponding to Ife = I(^ks) I{k,2) > ' ' ' , I{k,va) 



SNR1SNR2, 
by: 



■ , SNRr, 



SNR„ 



10Iog(%l)di? 



(36) 



In the update operator, because cr* is a function of the SNR, it becomes 



(4''^) in m- 



Note that the 9k in ([5]) should also be adjusted. Because weak signal infor- 
mation should be preserved by TBD algorithms, the Ik in ([SJ is the minimum 

of ife. 
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For the targets with unknown SNRs, it is assumed that the range of the 
SNRs of the targets is [SNRi, SNRh], corresponding to [I(k.i} , I{k,h)]- The I^^ 
could be normahzed for the region [Iik.i)^ I(k,h)]- The update operator and 
threshold are similar to those of multiple targets with different SNRs. 



6 Simulation 

6.1 Multiple targets miss distance 

The optimal sub-pattern assignment (OSPA) distance [15] between the esti- 
mated and true multitarget state is adopted here to estimate error. 

The OSPA distance is defined as follows [TS]. Let y) :— min(c, ||q — y||) 

for q, y G 3?^, and Up denote the set of permutations on {1,2, ... ,p} for any 
positive integer p. c is the cut-off value. Then, for Q = {qi,---,qjn} and 
Y = {yi, • • • ,y„}, 

J(=)(g,r) :=i ( min V d^^^ (qi, y.(i)) + c * (n - m) ] (37) 

if m < n; andJf'^) iQ,Y) -.^ d^"'^ [Y, Q) if m > n; and {Q,Y) = {¥, Q) = 
Q \i m = n — Q. 

For example, when m = 1 and n = 2, Qi = {gi} and Yi — {2/1,2/2}, the 
d^"'^ (Qi, Yi) = 0.5 * {min{min{\qi - yi\,c), min{\qi - 2/2|,c)) + c). 



6.2 Multiple targets with the same SNR 

Consider range cells in the interval [80000, 90000] with meters as the unit and 
Doppler cells in the interval [—400,-150] with meters/second as the unit. 
Nr — 200, Nd = 10 and A*";, = 1. The size of a range and Doppler cell is as 
follows: R=50 m and D=25 m/s. Lk = 2000 and Jk = 800 . The time interval 
T is equal to 1 s. 

At time step 1, the first target appears with the initial target state [ 89000 —200 
At step 10, a secondary target spawns from the first target, moving at -300 
m/s in the x-direction. The standard deviation of the noise is 0.25. 

In the first example, the SNRs of both targets are 9 dB. The measurements 
are shown in Fig. 3. Tracking with the shrinkage-PHD filter, we plot the PHD 
of the position and velocity versus time in Fig. 4 and 5. It is shown that 
the PHD of targets can be integrated over time by simultaneous detection 
and tracking. Observe that both trajectories are automatically initiated and 
tracked. In step 3, the first target can be tracked stably. At step 11, just 1 s 
after target spawning, the second target is already initiated and tracked. 

For the second example, we conducted simulations to compare the perfor- 
mance of the shrinkage-PHD filter (shrinkage-PHDF) the PHD filter (PHDF) 
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Fig. 4 PHD as a function of position and time step, using the shrink-PHD filter. 
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Fig. 5 PHD as a function of position and time step, using the shrink-PHD filter. 

[I6] , and the classical multitarget particle filtering method (MPF) [8] in track- 
ing accurately and consistently, and 50 Monte Carlo runs for the same scenario 
as used in the first example were performed. The states of targets are estimated 
by expectation-maximization (EM) algorithm. In addition, the maximum pos- 
sible number of targets is limited and known for MPF, but unlimited and 
unknown for the PHD filter and the shrinkage-PHDF. 
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To measure the estimation quantitatively, the estimation errors in terms of 
Monte Carlo averaged OSPA distance are shown in Figs. 6 and 7. For position 
estimation, the cut-off value c is five times the size of the range cell. For 
velocity estimation, the cut-off value c is two times the size of the velocity cell. 
It is shown that the shrinkage-PHDF estimates target position and velocity 
on all tracks with higher estimated precision than the MPF, because MPF 
requires a modeling setup to accommodate varying numbers of targets, which 
is very difficult to accomplish in reality.lt also shows that the performance of 
shrinkage-PHDF is better than PHDF, which is discussed further below. 
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Fig. 6 Comparison of the estimation errors of position given by the Monte Carlo-averaged 
OSPA miss distance with a measurement resolution R=50 m. 




Fig. 7 Comparison of the estimation errors of the velocity given by the Monte Carlo- 
averaged OSPA miss distance with a measurement resolution D=25 m/s. 
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To compare the shrinkage-PHDF and the PHDF clearly, Figs. 8 and 9 show 
the comparison of the mean OSPA for varying SNRs. The SNR of the targets 
is set to known values between 6 and 13 dB. When the value is between 6 and 
11 dB, the shrinkage-PHDF performs significantly better than the PHDF. 
When the value is large than 11 dB, they perform similarly. This finding 
illustrates the limitation of ordinary the PHDF for TBD: the PHD cannot 
sufficiently approximate the multitarget posterior probability when the SNR is 
low. The superiority of the shrinkage-PHDF takes advantage of the shrinkage 
operation is indicated, especially for low SNR scenarios. It is clear that for 
multitarget TBD, the shrinkage-PHDF performs better than the PHDF and 
MPF, especially when the SNR is low. 
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Fig. 8 Comparison of the estimation errors of position given by the time-averaged OSPA 
miss distance with a measurement resolution R=50 m for varying SNRs. 
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Fig. 9 Comparison of the estimation errors of the velocity given by the time-averaged OSPA 
miss distance with a measurement resolution D=25 m/s for varying SNRs 
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6.3 Multiple targets with different SNRs 

In this example, the targets with different SNRs (the first target is 8dB, the 
second target is 9dB and the last one is lOdB). At time stepl, the first target 
appears as in the first example. At step 7, the second target spawns form the 
first one, and the velocity in the x-direction is -300m/s. The third target new- 
birth far from the previous two targets at the time step 13 with the intital 
target state [89000 —250 O] .A comparison of the estimation errors of the 
position and the velocity of the three algorithms is shown in Fig. 10 and 11. 
In this scene, the targets vary not only the value of SNRs but also the average 
number of targets. Because the SNR of the first target is lower than which in 
section 6.2, the convergence velocity of the estimation first target is slower. 
After time step 7, comparing Figs. 6 and 10, Fig. 7 and 11, the varying SNRs 
and the number of the targets did not significantly influence the performance 
of the PHDF and Shrinkage-PHDF. In this example, the performance of MPF 
is better than that in the first example. The reason is the intensity of the 
targets is the most important characteristic to differentiate different targets 
for MPF method. In this way, the SNRs of all targets should be known for 
MPF method. 



-«-PHDF 

-e- Shrinkage-PHDF 
-A- MPF 

■ ■ ■ Cut-off value{5R) 




Fig. 10 Comparison of the estimation errors of the position given by the Monte Carlo- 
averaged OSPA miss distance for multiple targets with different SNRs (8dB, 9dB and lOdB). 



6.4 Multiple targets with unknown SNRs 



The PHD filter and the shrinkage-PHD filter are compared in Fig. 12 and 13. 
The sense is similar with that in the section 6.2, but the SNRs of the targets 
are unknown and the SNR of the second target changes from 9dB to lOdB. 
It is obvious that the performance of both algorithms is similar to what was 
presented in Figs. 6 and 7. After the 11th time step, the performance is even 
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Fig. 11 Comparison of the estimation errors of the velocity given by the Monte Carlo- 
averaged OSPA miss distance for multiple targets with different SNRs (8dB, 9dB and lOdB). 



better than which in the section 6.2, because the value of the SNR of the 
second target is increased. Therefore, whether the SNRs of muUiple targets 
are known has not significant influence on the performance of PHDF and 
Shrinkage-PHDF. In addition, MPF from 8 does not work when the value of 
the SNRs of the targets are unknown. 
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Fig. 12 Comparison of the estimation errors of the position given by the Monte Carlo- 
averaged OSPA miss distance for multiple targets with unknown SNRs (9dB and lOdB). 



6.5 A comparison of computation resources 

Table 3 shows the computation resources needed for the PHDF, the shrinkage- 
PHDF and MPF. All simulations were run on a PC with an Intel (R) Core 
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Fig. 13 Comparison of the estimation errors of the velocity given by the Monte Carlo- 
averaged OSPA miss distance for multiple targets with unknown SNRs (9dB and lOdB). 



Table 3 Computation resources 



Algorithm 


Shrinkage-PHDF 


PHDF 


MPF 


Computation(Second) 


28.66 


28.12 


84.20 



TM2 Duo CPU E4500 @ 2.20 GHz processor. The time listed is Monte Carlo- 
averaged, for a run in a scenario similar to that given in Section 6.2. The 
amount of calculation time required for the shrinkage-PHDF and the PHDF 
is much less than for MPF. The calculated amount time required for the 
shrinkage-PHDF is slightly greater than that required for the PHD filter. 



7 Conclusion 

In this paper, we used the RFS to model multitarget TBD measurements and 
to design an efficient Shrinkage-PHD filter for multitarget TBD. This filter 
accompanies with a shrinkage operation, and the optimal parameter for the 
shrinkage operation is obtained by an optimization procedure. Simulations 
show that the shrinkage-PHD filter takes very little time to detect new tar- 
gets, and it is sensitive to the variation in the number of targets. Moreover, 
our algorithm can track targets with high accuracy by taking advantage of 
the shrinkage operation. In addition, there is no restriction on the maximum 
number of possible targets or known value of SNRs of targets. In a word, by 
the application of Shrinkage-PHD filter, multiple targets could be tracked well 
in low SNR environments. 

In the future research, a challenge to be explored for the shrinkage-PHD 
filter will be determining how to model the measurement of extended targets 
in the framework of FISST. If this issue is resolved, it is believed that our new 
approach will be useful in future radar systems using the TBD technique. 
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